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Abstract 



The cosmic ray primary composition in the energy range between 10 and 10 
eV, i.e., around the "knee" of the primary spectrum, has been studied through the 
combined measurements of the EAS-TOP air shower array (2005 m a.s.L, 10 5 m 2 
collecting area) and the MACRO underground detector (963 m a.s.L, 3100 m w.e. of 
minimum rock overburden, 920 m 2 effective area) at the National Gran Sasso Lab- 
oratories. The used observables are the air shower size (N e ) measured by EAS-TOP 
and the muon number (N„) recorded by MACRO. The two detectors are separated 
on average by 1200 m of rock, and located at a respective zenith angle of about 
30°. The energy threshold at the surface for muons reaching the MACRO depth is 
approximately 1.3 TeV. Such muons are produced in the early stages of the shower 
development and in a kinematic region quite different from the one relevant for 
the usual — N e studies. The measurement leads to a primary composition be- 
coming heavier at the knee of the primary spectrum, the knee itself resulting from 
the steepening of the spectrum of a primary light component (p, He). The result 
confirms the ones reported from the observation of the low energy muons at the 
surface (typically in the GeV energy range), showing that the conclusions do not 
depend on the production region kinematics. Thus, the hadronic interaction model 
used (CORSIKA/QGSJET) provides consistent composition results from data re- 
lated to secondaries produced in a rapidity region exceeding the central one. Such 
an evolution of the composition in the knee region supports the "standard" galactic 
acceleration/propagation models that imply rigidity dependent breaks of the differ- 
ent components, and therefore breaks occurring at lower energies in the spectra of 
the light nuclei. 



1 Introduction 



The study of the primary cosmic ray composition and of its evolution with 
primary energy is the main tool in understanding the cosmic ray acceleration 
processes. In particular the energy range between 10 15 and 10 16 eV is char- 
acterized by breaks in the size spectra of the different Extensive Air Shower 
(EAS) components: electromagnetic (e.m.) [1], muon [2], Cherenkov light [3], 
and hadrons [4], which are therefore interpreted as a break in the primary 
energy spectrum. It is now recognized that the interpretation of such a fea- 
ture could provide a significant signature in understanding the galactic cosmic 
radiation [5], [6], [7], [8], [9]. 

Independent measurements based on the observation of the e.m. and GeV 
muon components [10], [11] lead to a composition becoming heavier in this 
energy region. The situation is more complex when other components are 
considered, thus showing that further information is needed from independent 
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observables (see, e.g., [12] and references therein). This is also useful to cross 
check the information, reduce the dependence on the hadron interaction model 
and particle propagation codes used, and to have better control of fluctuations 
in shower development, and therefore of event selection. 

At the National Gran Sasso Laboratories, we have developed a program of 
systematic study of the surface shower size measurements from EAS-TOP 
and the high energy muons [E 1 ^ = 1.3 TeV) measured deep underground 
(MACRO). Such muons originate from the decays of mesons produced in the 
first interactions of the incident primary in the atmosphere, and thus are from 
a quite different rapidity region than the GeV muons usually used for such 
analyses (xp > 0.1 or 0.2, the rapidity region being y — y^eam ~ —(4.5 — 5.0) 
at ^fs « 1000 TeV). The experiment provides therefore new data related to 
the first stages of the shower development, from secondaries produced beyond 
the central rapidity region. 

EAS-TOP and MACRO operated in coincidence in their respective final con- 
figurations for a live time of AT = 23,043 hours between November 25, 1992 
and May 8, 2000, corresponding to an exposure T ■ AT « 4 x 10 9 m 2 s sr. We 
present here an analysis of the full data set. Further details and partial results 
of the present work can be found in [13], [14] and [15]. 



2 The detectors 

The EAS-TOP array was located at Campo Imperatore (2005 m a.s.L, at 
about 30° from the vertical at the underground Gran Sasso Laboratories, cor- 
responding to an atmospheric depth of 930 g cm -2 ). Its e.m. detector (in which 
we are mainly interested in the present analysis) was built of 35 scintillator 
modules each 10 m 2 in area, resulting in a collecting area A « 10 5 m 2 . The ar- 
ray was fully efficient for N e > 10 5 . In the following analysis, we will use events 
with at least 7 neighboring detectors fired and a maximum particle density 
recorded by an inner module ("internal events"). The EAS-TOP reconstruc- 
tion capabilities of the EAS parameters for such events are: 10% above 
N e PS 10 5 for the shower size, and A9 0.9° for the arrival direction. The 
array and the reconstruction procedures are fully described in [16]. 

MACRO, in the underground Gran Sasso Laboratories at 963 m a.s.L, with 
3100 m w.e. of minimum rock overburden, was a large area multi-purpose 
apparatus designed to detect penetrating cosmic radiation. A detailed de- 
scription of the apparatus can be found in [17]. In this work we consider only 
muon tracks which have at least 4 aligned hits in both views of the horizon- 
tal streamer tube planes out of the 10 layers composing the lower half of the 
detector, which had dimensions 76.6 x 12 x 4.8 m 3 . The MACRO standard 
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reconstruction procedure [18] has been used, which provides an accuracy due 
to instrumental uncertainties and muon scattering in the rock of 0.95° for 
the muon arrival direction. The muon number is measured with an accuracy 
AiV M < 1 for multiplicities up to xs 10, and AN^ ss 1 for > 10; high 
multiplicity events have been scanned by eye to avoid possible misinterpreta- 
tions. 

The two experiments are separated by a rock thickness ranging from 1100 
to 1300 m, depending on the angle. The energy threshold at the surface for 
muons reaching the MACRO depth ranges from Ejj 1 = 1.3 TcV to E^ 1 = 1.8 
TeV within the effective area of EAS-TOP. 

The two experiments operated with independent triggering conditions, set to 
1) four nearby detectors fired for EAS-TOP (corresponding to a primary en- 
ergy threshold of about 100 TeV), and 2) a single muon in MACRO. Event co- 
incidence is established off-line, using the absolute time given by a GPS system 
with an accuracy better than 1 fis. The number of coincident events amounts 
to 28,160, of which 3,752 are EAS-TOP "internal events" (as defined above) 
and have shower size N e > 2 x 10 5 ; among them 409 have N e > 10 5,92 , i.e., are 
above the knee observed at the corresponding zenith angle [19]. We present 
here the analysis of such events, by using full simulations 1) of the detec- 
tors (based on GEANT [20]), 2) of the cascades in the atmosphere performed 
within the same framework as for the surface data (CORSIKA/QGSJET [21]), 
and 3) of the MUSIC code [22] for muon transport in the rock. Independent 
analyses from the two experiments separately are reported in [23], [19] and 
[11]. ' " ' ' 



3 Analysis and results 

3.1 The data 

The experimental quantities considered are the muon multiplicity distribu- 
tions (for > 1 as required by the coincidence trigger condition) in several 
intervals of shower sizes. We have chosen six intervals of shower sizes covering 
the region of the knee: 

5.20 < Log 10 (N e ) < 5.31 (1432 events), 5.31 < Log 10 (N e ) < 5.61 (2352 events), 
5.61 < Log 10 (N e ) < 5.92 (881 events), 5.92 < Log 10 (N e ) < 6.15 (252 events), 
6.15 < Log w (N e ) < 6.35 (106 events) and 6.35 < Log w (N e ) < 6.70 (42 
events). The experimental relative frequencies of the multiplicity distributions 
are shown in Fig.l. For further analysis, the data have been grouped in variable 
multiplicity bin sizes reported with their contents in Tables 1-6. 
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Table 1 



Relative frequency multiplicity distribution for size window: 5.20 < Logio(N e ) < 
5.31 (1432 events). The simulated distributions are also reported for the Light, 
Heavy components Monte Carlo (MC), and the resulting fit (see text). The number 
of digits is chosen in order to show the one event level. 
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Table 2 

As Table 1 for size window: 5.31 < Log w (N e ) < 5.61 (2352 events). 
3.2 The simulation 



We have simulated atmospheric showers in an energy range which includes 
the e.m. size values considered here (between 100 and 100,000 TeV/particle) 
and in an angular range exceeding the aperture of the coincidence experiment. 
Shower simulations have been performed with the QGSJET [24] hadronic in- 
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Table 3 

As Table 1 for size window: 5.61 < Log w (N e ) < 5.92 (881 events). 
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Table 4 -— ~ - "'" 
As Table 1 for size window: 5.92 < Log w (N e ) < 6.15 (252 events). 
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Table 5 

As Table 1 for size window: 6.15 < Logio(N e ) < 6.35 (106 events). 
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Table 6 

As Table 1 for size window: 6.35 < Log w (N e ) < 6.70 (42 events). 
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Fig. 1. Distributions of the relative frequencies of the detected underground muon 
multiplicities in the 6 selected size windows. See Tables 1-6. Notice the expected 
increasing relative frequency of high multiplicity events as a function of shower size. 



teraction model implemented in CORSIKA v5.62 [21] * . Primary particles 
have been sampled in a solid angle region of the order of the area encom- 
passing the surface array as seen from the underground detector. The solid 
angle corresponding to the selected angular window is Q = 0.0511 sr. All 
muons with energy > 1 TeV reaching the surface have been propagated 
through the rock down to the MACRO depth by means of the muon trans- 
port code MUSIC [22]; the accuracy of this transport code has been verified 
by comparing its results to those achieved with other Monte-Carlo simula- 
tions. Generated events having no muons surviving underground have been 
discarded, while those having at least one surviving muon have been folded 
with the underground detector simulation according to the following method, 
whose theoretical principles are discussed in [29]. We have considered an array 
of 39 (13 x 3) identical MACRO detectors adjacent to one another, covering 



* Simulations, using the other hadronic models in CORSIKA (DPMJET [25], 
HDPM [26] and SIBYLL [27]) have been performed, with reduced statistics, in order 
to verify the consistency of the procedure. The results are discussed in [15,28]. 
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an area of 230.7 x 158.2 m 2 . The shower axis is sampled over the horizontal 
area of the central MACRO, and all hit detectors are considered. For each 
hit detector, the full GMACRO (GEANT based) simulation of MACRO is 
invoked and is considered as a different event. For each of these events, when 
considered at the position of the "real" MACRO, the position of the shower 
core at the surface is recalculated, the particle densities on EAS-TOP counters 
are calculated and the trigger simulation is then invoked. Particle densities are 
obtained from the lateral distribution of the e.m. component of the shower as 
produced by CORSIKA (with the analytical "Nishimura-Kamata-Greisen" op- 
tion), taking into account the fluctuations of the number of particles hitting 
the detector modules and the full detectors' fluctuations [16]. If the trigger 
threshold is reached, the reconstructions of both EAS-TOP and MACRO are 
activated, thus producing results in the same format as the real data. The 
resulting events, after combining the simulated reconstructions of surface and 
underground detectors, are eventually stored as simulated coincidence events. 
Five samples with different nuclear masses have been generated: proton, He- 
lium, Nitrogen (CNO), Magnesium and Iron, all with the same spectral index 
7 = 2.62. Shower size bins have been chosen to be small enough so that no 
significant change in the shape of the muon multiplicity distributions in each 
bin is observed for different, extreme spectral indexes. A number of events 
exceeding the experimental statistics have been simulated in each size bin. 



3.3 The results 



The analysis is performed through independent fits of the experimental muon 
multiplicity distributions in the selected intervals of shower size. The simulated 
multiplicity distributions have been used as theoretical expectations for the 
individual components, and the relative weights are the fit parameters. 

The possibility that the experimental data could be reproduced with a single 
mass component can be easily excluded for the extreme (p or Fe) components, 
but also for medium mass primaries (e.g., A = 14): the obtained values of % 2 
are not satisfactory (too large, the point will be addressed at the end of the 
section). On the other hand, with a number of components larger than two we 
cannot achieve better solutions, since in all cases the minimization algorithm 
tends to force to zero the contribution of an intermediate component. This is 
mainly due to our limited statistics. For this reason we performed our analysis 
by considering only two components in the primary beam. We have tested 
two cases: a combination of p and Fe components, and a combination of two 
admixtures: a "Light" (L) and "Heavy" (H) one, built with equal fractions 
of p plus He and Mg plus Fe, respectively. Preliminary results from the p+Fe 
analysis have been presented in [15,28]. Here we describe the final analysis in 
terms of the L+H admixtures. 
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Table 7 

The fitted normalizations for the two components (L, H) as a function of size (notice 
that the two parameters are correlated, so that errors are not independent from one 
another). 

The fit has been performed in the six quoted size windows by minimizing the 
following expression for each multiplicity distribution: 



E 



(Nr - PlK l - p H N») 2 

°lexp + (PL^i,L) 2 + (PH^h) 



(1) 



where N^ xp is the number of events observed in the i th bin of multiplicity (with 
statistical uncertainty <Ti texp ), Nf 1 and are the numbers of simulated events 
in the same i th multiplicity bin from the L and H components, respectively, 
Pl and pu are the parameters (to be fitted) defining the fraction of each mass 
component contributing to the same multiplicity bin, and and a i H are 
the statistical errors of the simulation. Such an expression is close to that of a 
X 2 , although, in principle, it follows a different statistics, and in the following 
we shall refer to it as if it were a genuine x 2 - The values of the parameters p^ 
and ph obtained from the minimizations are given in Table 7. The progressive 
decrease of the "Light" component in favor of the "Heavy" one is visible and 
significant at the level of 2 standard deviations: the average pl value is 0.70 
± 0.04 below the observed knee in size (Logio(iV e ) = 5.92), and 0.28 ± 0.17 
above. By normalizing p^ and pu to the observed number of coincident events 
in each size bin (see Tables 1-6) we obtain the contribution to the measured 
size spectrum of each component. 

In Fig. 2 the multiplicity distributions are shown for the four most relevant 
size windows, together with the expected L and H components, and their best 
fit combination. 

Regarding the shapes of the multiplicity distributions, it is interesting to re- 
mark that they cannot be described by simple single laws, and show some 
structure; this is evident in the data and in the simulated Heavy components, 
however less so in the simulated Light ones. The origin of such structure is 
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Fig. 2. Distributions of the relative frequencies of the detected underground muon 
multiplicities (black points) together with the predictions for L (open triangles) and 
H (open stars) admixtures in the QGSJET interaction model, and the (L+H) fit 
(open squares). 



entirely geometric and due to the interplay between the typical size of muon 
bundles with the two length scales of the MACRO detector. Small bundle sizes 
can be entirely contained in the detector while, when the size increases, this 
becomes impossible along the width of the detector. Bundles of even larger size 
exceed also the length of MACRO. This fact is well taken into account by the 
simulation, and in fact the fit reproduces correctly this change of structure, 
which is typical of large bundles (i.e., high energies and large masses). The 
effect is evident when comparing with a single component fit, say the CNO 
group that has an intermediate average atomic number. The results of the fit 
are presented in Fig. 3. CNO primaries alone provide good fits in the higher 
size bins (due to the limited statistics), but below and just above the knee 
at Log w (N e ) = 5.92, the large x 2 values indicate the failure to reproduce the 
shape of the multiplicity distribution (see Table 8). 
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Table 8 

The x 2 values resulting from the fits to the CNO (A=14) component alone, as a 
function of shower size. 
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Fig. 3. Distributions of the relative frequencies of the detected underground muon 
multiplicities (black symbols) together with the results of the fits for the QGSJET 
interaction model (L+H) (open squares), compared with a fit with the CNO com- 
ponent only (open triangles). 

3.4 Interpretation of the data 



For a given size window, the contribution of each primary mass group derives 
from a different energy region: the higher the mass number, the higher the 
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corresponding energy. The size-energy mapping far from shower maximum is 
model dependent, and in our analysis is based on CORSIKA/QGSJET. From 
the full simulation chain we also calculate the probabilities e a (E,AiN e ) for a 
primary belonging to mass group a (a = L, H) and of energy E to give a 
coincident event in the i th size window AiN e . To evaluate the average mass 
composition we use a logarithmic energy binning (3 bins per energy decade), 
starting from 100 TeV/nucleus. From the simulation we obtain the number 
of events (n"(AjiV e )) that a primary of mass group a will produce in the j th 
energy bin, when the detected size is in the windows AjiV e . Therefore the total 
number of events that the primary mass group a produces in the size window 
Aj7V e is the sum of n"(A;iV e ) over the energy bins. 

We require that the number of experimentally observed events in the size 
window AiiV e be equal to: 

N e ^(A t N e )=p L (A l N e )J2nf(A l N e )+p H (A l N e )J2nf(A t N e ) (2) 

3 3 

where p L and p H are the fit coefficients for the given size window AjiV e . These 
are normalized, so that p L = 1 — p H in each size window. This leaves an 
overall renormalization factor K free in order to satisfy eq.2, so we obtain the 
renormalized quantities n* a = K n". The corrected estimated number of 
primaries of mass group a for each size window belonging to energy bin j can 
thus be obtained by applying the efficiencies e a (Ej, AjiV e ): 

m?(A,Ag = p a (A i iV e )nf (A i JV e )/e Q (^-, A,N e ) (3) 

Then, since the j th energy bin may receive contributions from different size 
windows, we have to sum over % (the size window index): 

Mf = 5>J(A,iV e ) = ^(^Mn^iAM/e^E^AM (4) 

i i 

Mf and Mf provide estimates of the energy spectra of the L and H mass 
groups, presented in Fig. 4. There we plot the spectra starting from 10 3 TeV 
since with our selection of size, this is the energy at which the heaviest com- 
ponent has reached a significant triggering efficiency 

A steepening by about A7 = 0.7 ± 0.4 of the light mass group spectrum 
just at the knee (~4xl0 15 eV) is observed, assuming power law behaviors 
crossing at the knee position. Although these distributions cannot be used 
to obtain a direct representation of the actual cosmic ray spectrum, due to 
the two mass groups schematization and the choices of their components, the 
relative proportion of "Light" and "Heavy" admixtures turns out to be quite 
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Fig. 4. Energy spectra estimates for L and H admixtures. 



stable with respect to the mentioned parameters; within this approximation, 
the resulting all-particle spectrum would show a A7 = 0.4 ± 0.1. 

We make use of the values of Mf and Mf so obtained to compute the mean 
value of the natural logarithm of the primary mass (< In A >) as a function 
of energy: 



< \nA{Ej) > 



In A L Mf + In A H Mf 

Mf + Mf 



(5) 



with lnA L = 0.5(lnA p + lnA He ), and lnA H = 0.5(lnA M v + lnA Fe ). The uncer- 
tainty on < In A(Ej) > has been obtained by propagating the uncertainties 
on the fit coefficients. The result is reported in Fig. 5 together with the results 
of KASCADE [10] and EAS-TOP alone [11], where these analyses has been 
performed using e.m. size and GeV muons detected at surface level. The good 
agreement shows that the results do not depend on the selected muon energy. 
The < In A(Ej) > obtained by MACRO alone [23], on the basis of the HEMAS 
Monte Carlo code [30], has a milder energy dependence and appears to be in 
contrast with those presented here above Logio(E) > 4.2. In our opinion this is 
due to a weakness of the HEMAS model, based on parameterizations of UA5 
results [31]. The possible shortcomings of the HEMAS model were already 
discussed in [32,33]. 

The allowed region for < InA(E) > obtained from our analysis is also consis- 
tent with the theoretical expectations from refs. [9], [34], and [35]. 
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Fig. 5. < In A > vs primary energy (continuous line). The hatched areas represent 
the 68% uncertainty range due to the statistical error. We also superimpose the 
results of KASCADE [10] and of EAS-TOP [11] (open squares). 



4 Conclusions 



The analysis N e -N^ ey events collected by the MACRO/EAS-TOP Collabora- 
tion at the Gran Sasso Laboratories points to a primary composition becoming 
heavier around the knee of the primary spectrum (i.e., in the energy region 
10 15 — 10 16 eV). The result is in good agreement with the measurements of 
other experiments based on the observation of the e.m. and muon compo- 
nents at ground level. The muon energies detected in the present experiment 
are however about three orders of magnitude larger than in previous N e -N M 
experiments, and therefore the parent pions are produced in a different kine- 
matic region (at the edges of the fragmentation region, rather than in the 
central one) and in the first stages of the cascade development. A good overall 
consistency of the interaction model used (CORSIKA/QGSJET) in describ- 
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ing the yield of secondaries over a wide rapidity region is thus obtained . 
The present data explain therefore the observed knee in the cosmic ray pri- 
mary spectrum as due to the steepening of the spectrum of a light component 
(p, He) at E « 4 x 10 15 eV, of A7 = 0.7 ± 0.4. Such an effect can be inter- 
preted in the "standard" framework of the acceleration/propagation processes 
of galactic cosmic radiation that predict, as a general feature, rigidity depen- 
dent breaks for the different nuclei, and therefore appearing at lower energies 
for the lighter ones. 
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